function photon = MonteCarlo(mu_a, mu_s, n1, n2, n, g)

mu_t=mu_a + mu_s;

Wbdy=0;
photon=[];
T=[];
photon_xyz = [0 0 0];    % Initial position of the photon.
photon_dir = [0 0 1];    % Initial direction of the photon.
S = 0;                   % Initial step size of photon.
W = 1;                   % Weight of photon.
ksi = 0;
theta_new = 0;
phi_new = 0;
Rsp=(n1-n2)^2/(n1+n2)^2;    %Specular reflectance in normal incidence
W=W*(1-Rsp);                %Transmitted intensity at surface
 
 % Loop until photon dies, reaches a boundary, or is transmitted.
 while 1
   
    %Store the location,direction and energy of the photon
    photon=[photon;photon_xyz,photon_dir,W];   
    
    % ksi value (assign random values from a uniform distrobution ranging
    % from 0->1)
    ksi1=rand();
    ksi2=rand();
    ksi3=rand();
    
    % Find the step size for this iteration.
    S = -log(ksi1);
    
    % Update the location of the photon based on it's scattering.
    % (e.g. x' = x + step_size*direction)
    temp=photon_xyz;
    photon_xyz(1,1) = photon_xyz(1,1) + S*photon_dir(1,1);
    photon_xyz(1,2) = photon_xyz(1,2) + S*photon_dir(1,2);
    photon_xyz(1,3) = photon_xyz(1,3) + S*photon_dir(1,3);
    W = W - mu_a/mu_t*W;
        
    % Get the new theta angle.  
    theta_new = theta(g, ksi2);
    
    % Get the new phi angle.
    phi_new = phi(ksi3);
    
    % Get the new direction cosines.
    photon_dir = mu(photon_dir, theta_new, phi_new);
    
    % Check to see if photon is inside the medium still
    % by calling function 'check_medium(...)'.  Returns
    % true (1) if photon is still within the medium.
    if((check_medium(photon_xyz))==1)
    
        % Check to see if photon has struck the boundary.  If it has
        % we absorb the energy and kill the photon.
        if((check_bounds(photon_xyz))==0)
            % FIX ME:
            %   Added this, should it be here?  It seems if it makes it to the
            %   boundary we should record where that happens and the
            %   energy it had when that happened.
            photon=[photon;photon_xyz,photon_dir,W]; 
            
            %Wbdy=W; %Absorb energy of photon at boundary
            %display('photon absorbed in a boundary');
            break;
        end %Do nothing and go ahead if photon didnt strike any boundary
    
            % Perform roulette to see if photon regains energy based on
            % probability or if the photon is killed off.
            if (W < 0.001)  %Is the energy is below the given threshold?
                if (1/10 >= rand()) %Did the photon survive roulette?
                    W = 10 * W;
                    %FIX ME:
                    %   Added this, shouldn't we just continue the loop from this
                    %   point?
                    continue;
                   
                    %display('Photon ragained energy from roulette');
                else %If photon didnt survive roulette
                    % FIX ME:
                    %   Added this. If the photon doesn't survive
                    %   roulette should we record where it died off?  Also
                    %   its energy when that happened?
                    photon=[photon;photon_xyz,photon_dir,W]; 
                    
                    W = 0;  %Kill the photon by setting its energy to zero.
                    %display('Photon terminated at roulette');
                    break;  %If the photon was terminated we break from the loop.
                end %End of roulette
            end % End of energy threshold check (Do nothing and carry on as usual if energy is above threshold)
      else %Otherwise the photon is outside the medium
                photon_xyz=temp;
                [R, alpha_t] = refl(photon_xyz, photon_dir, n1, n2);
                ksi4=rand();
                if(ksi4<=R)%If Total Internal Reflectance probabilistically happens update the direction.
                    photon_dir(1,3)=-photon_dir(1,3);
                    %display('Photon totally internally reflected');
                else %TIR doesnt happen (i.e. the photon is transmitted)
                    T=[T; alpha_t W];
                    %display('photon transmitted');
                    break;
                end %End of reflectance check
       end%End of photon location in medium check
        
 end %End loop
    
    
    
    
    
    
    

    
    
    